STA 360 Lab 4: Introduction to Hamiltonian Monte Carlo

Overview

By now you have encountered at least two procedures for generating samples from a distribution over model parameters conditional on observed data (the posterior distribution of your model parameters). The first method, direct sampling, assumes that we can explicitly write down the form of the posterior density and/or distribution function, which means that we can generate samples from it with calls to standard R functions. The second method, Gibbs sampling, generates samples from a posterior distribution by building a Markov Chain whose stationary distribution is the desired posterior.

The benefit of Gibbs sampling is that it allows us to sample from distributions whose form we may not be able to write down. The drawback is that, unlike samples directly from the posterior, samples from the Gibbs sampler eventually behave as if they were drawn from the posterior distribution: now we have to be concerned about exactly how and when this convergence occurs. Later in the course you will learn about diagnostics for Gibbs samplers and some methods for monitoring convergence of Markov chains. You will write your own Gibbs samplers and implement some of these methods. You will also learn that the Gibbs sampler is just a special case of a class of procedures for sampling from a posterior distribution, which fall under the title Markov Chain Monte Carlo (MCMC).

Before you see those details, we will spend time in this lab building intuition for the mechanics of Markov Chain Monte Carlo methods and we will use that intuition to understand the advanced methods Stan is running under the hood to generate its samples.

MCMC and Gibbs sampling

One of the key intuitions to build when thinking about MCMC methods is that of “space exploration.” In the case of two-parameter models parametrized by \(\theta = [\theta_1, \theta_2] \in \mathbb{R}^2\), we can think of MCMC methods as sampling from the 2-d real-valued plane according to our target density

\[ p(\theta_1, \theta_2|X_1, \dots, X_n) \]

MCMC methods like the Gibbs sampler will move iteratively through the plane, transitioning from point to point and producing a set of \(S\) samples

\[ [\theta^{(1)}_1, \theta^{(1)}_2], \dots, [\theta^{(S)}_1, \theta^{(S)}_2] \]

When the target density is nicely behaved, the exploration of a Gibbs sampler is also “nice.” Regions of high density are visited with higher frequency than regions of lower density. All regions of the plane are both in theory and in practice accessible with only a few transition steps. An example of a nice density on \(\mathbb{R}^2\) is a bivariate normal with a modest degree of dependence between its two components \(\theta_1, \theta_2\):

What’s an example of a “not nice” density for moving through the plane with Gibbs transition steps? Let’s look at an example of a bivariate density that would likely give the Gibbs sampler some trouble:

The problem here is that there are regions with very high density (the red peaks in the middle) right next to regions of very low density. It will take many transition steps for the Markov Chain to reach the peaks, and will also take many transitions to get out of the neighborhood of the peaked regions once it moves near it.

While the asymptotic properties of the Markov chain still hold, the sticking behavior makes it difficult to get accurate posterior summaries without taking an infinite number of samples and waiting an infinitely long time for them. As with many other problems in statistics, this problem often gets worse as the number of dimensions in the parameter space increases.

Hamiltonian Monte Carlo (HMC)

Even less extreme examples than the one shown above will cause Gibbs samplers to explore the sample space more slowly than we’d like. Let’s look at a simpler case and write the code for our own Gibbs sampler. Later we will compare it to Stan’s behavior on the same problem and discuss the benefits of using HMC.

Consider data generated from a bivariate normal distribution with mean parameters \(\theta_1, \theta_2\) and known covariance matrix \(\Sigma\) and suppose we place independent normal priors on \(\theta_1, \theta_2\):

\[ \begin{aligned} X_1, \dots, X_n &\sim N \left(\theta, \Sigma \right) \\ \theta_j &\sim N(0, 1)~~~~~~j=1,2 \end{aligned} \]

Suppose too that the covariance matrix \(\Sigma\) has a specific form with unit-variance marginal distributions and known correlation parameter \(\rho\):

\[ \Sigma = \left[\begin{array}{cc} 1 & \rho \\ \rho & 1 \end{array}\right] \]

In this course you have derived (or will derive) the multivariate posterior density \(p(\theta | \Sigma, X)\). However, for the purposes of this exercise, we will sample from the full conditional densities

\[ p(\theta_1 | X_1, \dots, X_n, \rho, \theta_2) \\ p(\theta_2 | X_1, \dots, X_n, \rho, \theta_1) \]

Before building the Gibbs sampler to make inferences on \(\theta_1, \theta_2\), first answer these questions:

  1. What is the conditional density \(p(\theta_1 | X_1, \dots, X_n, \rho, \theta_2)\)?
  2. What is the conditional density \(p(\theta_2 | X_1, \dots, X_n, \rho, \theta_1)\)?

Now that you have derived the full condition densities, write some code to implement a Gibbs sampler for this model.

With the Gibbs sampling code in hand, let’s generate samples from the posterior distribution of \(\theta_1, \theta_2\) with \(\rho = 0.2\). We’ll do the same using Stan.

The Gibbs sampling results and the HMC results look pretty similar! What happens when \(\rho = 0.995\)?

Please answer these questions:

  1. How do the results of the Gibbs sampler differ from those obtained from HMC?
  2. Why do the samples from the Gibbs sampler exhibit this behavior?

STA 360: Bayesian Inference and Modern Statistical Methods

27 July, 2019